
** report monthly n, test rate, positivity, by group

local glist no_icli cancer labor_delivery ami stroke fracture wound ///
	appendicitis vehicle_accident other_accident
	
qui foreach group of local glist  {
	
	
	* load the data
	if "`group'" ~= "no_icli" use health/inp_clean if `group'==1, clear
	if "`group'" == "no_icli" use health/inp_clean if any_icli==0, clear
	keep if admit_date >= mdy(3,13,2020)
	gen month = month(admit_date)
	
	collapse (sum) test_tight pos_tight (count) n=test_tight, by(month)
	
	gen group = "`group'"
	save health/rates_`group', replace 
	
	
	gen test_rate = string(test_tight / n, "%4.3f") + " & "
	gen conf_pos_rate = string(pos_tight / n, "%4.3f") + " & "
	gen pos_rate = string(pos_tight / test_tight , "%4.3f") + "\\ "

	keep if month >= 3
	gen 		months = "March &" if month == 3
	replace months = "April &" if month == 4
	replace months = "May &" if month == 5
	replace months = "June &" if month == 6
	replace months = "July &" if month == 7
	replace months = "August &" if month == 8
	replace months = "September &" if month == 9
	replace months = "October &" if month == 10
	replace months = "November &" if month == 11
	replace months = "December &" if month == 12
	
	gen ns = string(n, "%10.0fc") + " & " 
	noi di ""
	noi di "`group'"
	noi list months ns test_rate conf_pos_rate pos_rate , noobs clean
	outsheet months ns test_rate conf_pos_rate pos_rate ///
		using tables/bounds_`group'.tex, replace noquote nonames delimit(" ")
	
}

********************************************************************************
********************************************************************************
** aggregate across months
********************************************************************************
********************************************************************************
local glist cancer labor_delivery ami stroke fracture wound ///
	appendicitis vehicle_accident other_accident
clear
foreach group of local glist {
	append using health/rates_`group'
}

keep if month >= 3 
desc

collapse (sum) test_tight pos_tight n, by(group) 
gen test_rate = test_tight/n
gen pos_rate = pos_tight/test_tight
gen lb = pos_tight/n

list group n test_rate lb pos_rate

** calculate se of proportions
	gen test_rate_se = sqrt(test_rate*(1-test_rate)/n)
	gen pos_rate_se = sqrt(pos_rate*(1-pos_rate)/test_tight)
	gen lb_se = sqrt(lb*(1-lb)/n)

** make nicely formatted table with confidence bands 
	rename (test_rate n) (test_rate1 n1)
	rename lb lb1
	gen lb2 = lb1 - 1.96*lb_se
	rename pos_rate ub1
	gen ub2 = ub1 + 1.96*pos_rate_se
	
	reshape long n test_rate lb ub, i(group) j(rr)

	
	** format strings for numbers
	tostring test_rate lb ub, replace force
	gen N = string(n, "%10.0fc") 
	drop n
	rename N n
	replace n = "" if rr == 2
	
	** form confidence intervals
	foreach v in test_rate lb ub {
		replace `v' = substr(`v', 1,4)
	}
	gen 		bound = "[" + lb + ", " + ub + "] \\" if rr==1
	replace bound = "(" + lb + ", " + ub + ") \\" if rr==2
	
	** blank out second row for name, n, test rate
	foreach v in group n test_rate {
		replace `v' = "" if rr == 2
	} 
	
	** rename
	replace group = "AMI" if group == "ami"
	replace group = "Appendicitis" if group == "appendicitis"
	replace group = "Cancer" if group == "cancer"
	replace group = "Fracture" if group == "fracture"
	replace group = "Labor/delivery" if group == "labor_delivery"
	replace group = "Other accident" if group == "other_accident"
	replace group = "Stroke" if group == "stroke"
	replace group = "Vehicle accident" if group == "vehicle_accident"
	replace group = "Wound" if group == "wound"
	
	** delimit 
	foreach v in group n test_rate {
		replace `v' = `v' + "  & " 
	}

	list group n test_rate bound , noobs clean
	outsheet group n test_rate bound using tables/groups_pool_time.tex, ///
		replace delimit(" ") nonames noquote





********************************************************************************
********************************************************************************
** Make figures
********************************************************************************
********************************************************************************
local glist cancer labor_delivery ami stroke fracture wound ///
	appendicitis vehicle_accident other_accident
clear
foreach group of local glist {
	append using health/rates_`group'
}

keep if month >= 3 
describe

replace group = "ami_stroke" if group == "ami" | group == "stroke"
replace group = "injury" if group == "wound" | group == "fracture"
replace group = "app_vehicle" if group == "appendicitis" | group == "vehicle_accident"

collapse (sum) test_tight pos_tight n, by(group month)

gen test_rate = test_tight/n
gen pos_rate = pos_tight / test_tight
gen conf_pos_rate = pos_tight/n

gen label = "Injury" if group == "injury" & month == 12
replace label = "App./vehicle" if group == "app_vehicle" & month == 12
replace label = "AMI/Stroke" if group == "ami_stroke" & month == 12
replace label = "Oth. accidents" if group == "other_accident" & month == 12
replace label = "Labor/delivery" if group == "labor_delivery" & month == 12
replace label = "Cancer" if group == "cancer" & month == 12

# delimit ;
twoway

	(connected test_rate month if group == "injury", color(black) lwidth(thick)
		mlabel(label) mlabpos(2) mlabcolor(black) mlabsize(medlarge))	
	(connected test_rate month if group == "app_vehicle", color(gs5) lwidth(thick) 
		mlabel(label) mlabpos(4) mlabcolor(gs5)	 lpattern(dash) mlabsize(medlarge))
	(connected test_rate month if group == "ami_stroke", color(gs9) lwidth(thick) 
		mlabel(label) mlabpos(4) mlabcolor(gs9) mlabsize(medlarge) )	
	(connected test_rate month if group == "other_accident", color(gs12) lwidth(thick) 
		mlabel(label) mlabpos(2) mlabcolor(gs12) mlabsize(medlarge))	
	
	(connected test_rate month if group == "labor_delivery", color(black) lwidth(thick)
		mlabel(label) mlabpos(2) mlabcolor(black)  lpattern(dash) mlabsize(medlarge))
	(connected test_rate month if group == "cancer", color(gs9) lwidth(thick) 
		mlabel(label) mlabpos(4) mlabcolor(gs9)	 mlabsize(medlarge) )
	,
	legend(off)
	xsize(7) ysize(5)
	graphregion(color(white) margin(r+25))
	xlabel(3 "Mar" 4 "Apr" 5 "May" 6 "Jun"
		7 "Jul" 8 "Aug" 9 "Sep" 10 "Oct" 11 "Nov" 12 "Dec")
	ytitle("") xtitle("")
;
# delimit cr
graph export figures/test_rate_narrow_group.pdf, replace 

** spike plot showing bound by monnth and group
gen 		gn = 1 if group == "cancer"
replace gn = 2 if group == "app_vehicle"
replace gn = 3 if group == "labor_delivery"
replace gn = 4 if group == "injury"
replace gn = 5 if group == "other_accident"
replace gn = 6 if group == "ami_stroke"

forvalues month = 3/11 {
	
	local max .18
	local oo ""
	if `month' == 3 {
		local mname "March"
		local oo "(off scale)"
		local max .3 
	}
	if `month' == 4 local mname "April"
	if `month' == 5 local mname "May"
	if `month' == 6 local mname "June"
	if `month' == 7 local mname "July"
	if `month' == 8 local mname "August"
	if `month' == 9 local mname "September"
	if `month' == 10 local mname "October"
	if `month' == 11 local mname "November"
	if `month' == 12 local mname "December"
	
	
	# delimit ;
	twoway
		(pccapsym pos_rate gn conf_pos_rate gn if month == `month')
		,
		xlabel(1 "Cancer" 2 "App./vehicle" 3 "Labor/delivery" 4 "Injury" 
			5 "Oth. accidents"  6 "AMI/Stroke", angle(45) labsize(medlarge))
		graphregion(color(white))
		xtitle("")
		ytitle("")
		text(`max' 6 "`mname'" "`oo'", place(w) size(large) justification(right))
		name(g`month', replace)
		ylabel(0(.06)`max')
		xsize(2) ysize(3)
	;
	# delimit cr
}


graph combine g3 g4 g5 g6 g7 g8 g9 g10 g11 , ///
	cols(3) xsize(6) ysize(9) graphregion(color(white))
graph export figures/bounds_by_narrow_group.pdf, replace 

sort month conf_pos_rate
list month group n test_rate conf_pos_rate pos_rate , sep(6)

window manage close graph _all
